Finite difference methods for approximating Heaviside functions

We present a finite difference method for discretizing a Heaviside function H(u(~x)), where u is a level set function u : Rn 7→ R that is positive on a bounded region Ω ⊂ R. There are two variants of our algorithm, both of which are adapted from finite difference methods that we proposed for discretizing delta functions in [13–15]. We consider our approximate Heaviside functions as they are used to approximate integrals over Ω. We prove that our first approximate Heaviside function leads to second order accurate quadrature algorithms. Numerical experiments verify this second order accuracy. For our second algorithm, numerical experiments indicate at least third order accuracy if the integrand f and ∂Ω are sufficiently smooth. Numerical experiments also indicate that our approximations are effective when used to discretize certain singular source terms in partial differential equations. We mostly focus on smooth f and u. By this we mean that f should be smooth in a neighborhood of Ω, and u should be smooth in a neighborhood of ∂Ω. However, our algorithms still give reasonable results if either f or u has jumps in its derivatives. Numerical experiments indicate approximately second order accuracy for both algorithms if the regularity of the data is reduced in this way. Numerical experiments indicate that dependence on the placement of Ω with respect to the grid is quite small for our algorithms. Specifically, a grid shift results in an O(hp) change in the computed solution, where p is the observed rate of convergence.